##########################################################################################
## https://sunduanchen.github.io/Scissor/vignettes/Scissor_Tutorial.html
library(data.table)
library(optparse)
library(ggplot2)
library(dplyr)
library(patchwork)
library(Seurat)

##########################################################################################

option_list <- list(
    make_option(c("--single_cell_file"), type = "character"),
    make_option(c("--gene"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    gene <- "MUC5AC"
    single_cell_file <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6/Scissor_STAD_MUC6_mutation.IM.CellRate.all.RData")
    out_path <- paste0("~/20220915_gastric_multiple/dna_combinePublic/images/singleCell_MUC6_vln/")

}

##########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

out_path <- opt$out_path
gene <- opt$gene
single_cell_file <- opt$single_cell_file

dir.create(out_path , recursive = T)

##########################################################################################

sc_dataset_all <- load(single_cell_file, verbose = F)
sc_dataset$celltype[which(sc_dataset$celltype=="Pit" & sc_dataset$scissor == 1) ] <- "Pit_Mut"
## Pit其它细胞
sc_dataset$celltype[which(sc_dataset$celltype=="Pit" & sc_dataset$scissor != 1) ] <- "Pit_Other"

sc_dataset_all <- sc_dataset
##Idents函数定义要取的类别
Idents(sc_dataset_all) <- "orig.ident"
sc_dataset_all <- subset(sc_dataset_all , idents=c("JZ732P"))

##########################################################################################
cell_order <- c("Enterocytes" , "Pit_Mut" , "Pit_Other" , "Endocrine" , "Goblet" , "Neck" , "Chief" , "Parietal")

Idents(sc_dataset_all) <- "celltype" 
## 细胞比例不看tumor细胞的太少了
plot <- VlnPlot(
    subset( sc_dataset_all , idents=cell_order ), 
    features = c(gene),
    pt.size = 0)

out_name <- paste0( out_path , "/" , gene , ".vln.IM.pdf"  )
ggsave(file=out_name,plot=plot,width=8,height=5)

##########################################################################################
## 肠化中MUC6重新画图

sc_dataset_all <- FindNeighbors(sc_dataset_all, dims = 1:10)
sc_dataset_all <- RunUMAP(object = sc_dataset_all, dims = 1:10)

plot <- FeaturePlot(object=sc_dataset_all,features=c("MUC6"),cols=c('grey','red'),label=F,pt.size=0.3)
out_name <- paste0( out_path , "/MUC6.FeaturePlot.IM.pdf"  )
ggsave(file=out_name,plot=plot,width=5,height=5)

plot <- FeaturePlot(object=sc_dataset_all,features=c("GKN1", "GKN2"),cols=c('grey','red'),label=F,pt.size=0.3)
out_name <- paste0( out_path , "/GKN1_GKN2.FeaturePlot.IM.pdf"  )
ggsave(file=out_name,plot=plot,width=10,height=5)

plot <- FeaturePlot(object=sc_dataset_all,features=c(gene),cols=c('grey','red'),label=F,pt.size=0.3)
out_name <- paste0( out_path , "/" , gene , ".FeaturePlot.IM.pdf"  )
ggsave(file=out_name,plot=plot,width=5,height=5)

